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Abstract 

This paper addresses the question of how population diffusion affects the formation of the spatial 
patterns in the spatial epidemic model by Turing mechanisms. In particular, we present theoretical 
analysis to results of the numerical simulations in two dimensions. Moreover, there is a critical 
value for the system within the linear regime. Below the critical value the spatial patterns are 
impermanent, whereas above it stationary spot and stripe patterns can coexist over time. We have 
observed the striking formation of spatial patterns during the evolution, but the isolated ordered 
spot patterns don't emerge in the space. 



I. INTRODUCTION 



The dynamics of spontaneous spatial pattern formation, first introduced to biology by 
Turing jl] five decades ago, has recently been attracting attention in many subfields of 



biology to describe various phenomena. Non-equilibrium labyrin 
in chemical reaction-diffusion systems with a Turing instability 



diffusion systems 
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;hine patterns are observed 
2J and in bistable reaction- 



Such dynamic patterns in a two dimensional space have recently 



been introduced into ecology 



3, H, Q, H 



]. In the past few years, geophysical patterns over 



a wide range of scales for the vegetation have been presented and studied in the Refs. 
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id, 



14j using the Turing mechanisms. 



In the epidemiology, one of the central goals of mathematical epidemiology is to pre- 
dict in populations how diseases transmit in the space. For instance, the SARS epidemic 
spreads through 12 countries within a few weeks. The classical epidemic SIR model de- 
scribes the infection and recovery process in terms of three ordinary differential equations 
for susceptibles (S), infected (I), and recovered (R), which has been studied by many re- 
searchers 15, If], [l?], [3] and the reference cited therein. These systems depend mainly on 
two parameters, the infection rate and the recovery rate. 



A growing body of work reports on the ro 



in the host population structure 
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e of spatial patterns on evolutionary processes 
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25]. Recent studies have shown 
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28|]. More 



large-scale spatiotemporal patterns in measles [26| and dengue fever (DF) 
dramatically the wave is often caused by the diffusion (or invasion) of virus within the 
populations in a given spatial region, thus generating periodic infection, which has been 
observed in the occurrence of dengue hemorrhagic fever (DHF) in Thailand 2j|. Existing 
theoretical work on pathogen evolution and spatial pattern formation has focused on a 
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25]. 



model in which local invasion to the susceptible hosts plays a central role [20l . 
Projections of the spatial spread of an epidemic and the interactions of human movement at 
multiple levels with a response protocol will facilitate the assessment of policy alternatives. 



3Q|, 



3l|. 



Spatially-explicit models are necessary to evaluate the efficacy of movement controls 
A wide variety of methods have been used for the study of spatially structured epidemics, 
such as cellular automata 
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equations 
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34j, networks 
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361 ] . metapopulations [37J, 138], diffusion 



41] , and integro-differential equations, which are useful tools in the study 



of geographic epidemic spread. In particular, spatial models can be used to estimate the 



2 



formation of spatial patterns in large-scale and the transmission velocity of diseases, and in 
turn guide policy decisions. 

This paper addresses how diffusive contacts and diffusive movement affect the formation 
of spatial patterns in two dimensions. The diffusion term is from the earlier work that 
tracing back to Fisher and Kolmogorov. Noble applied diffusion theory to the spread of 
bubonic plague in Europe Q|. Noble's model relies on the assumptions that disease trans- 
mits through interactions between dispersing individuals and infected individuals move in 
uncorrelated random walks. In light of the Turing theoretical and study of recent spatial 
models, we investigate the formation of spatial patterns in the spatial SIR model based on 



the study of non-spatial SIR model with constant removal rate of the infectives [15]. 



A — dS— XSI, (la) 



II. MODEL 

A. Basic model 

We consider, as the basic model, the following Susceptible- Infected-Recovery (SIR) model 

dS 

^ = XSI — (d + 7)1 - h(I), (lb) 

df? 

— = 7 J + h (I)- dR, (lc) 

where S(t), I(t), and R(t) denote the numbers of susceptible, infective, and recovered indi- 
viduals at time t, respectively. A is the recruitment rate of the population (such as growth 
rate of average population size, a recover becomes an susceptible, immigrant and so on), d 
is the natural death rate of the population, 7 is the natural recovery rate of the infective 
individuals, and A is a measure of the transmission efficiency of the disease from suscepti- 
bles to infectives. In Eq. p]) , h(I) is the removal rate of infective individuals due to the 
treatment. We suppose that the treated infectives become recovered when they are treated 
in treatment sites. We also suppose that 

r, for / > 0, 
h(I) = { (2) 

0, for 1 = 0, 
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where r > is constant and represents the capacity of treatment for infectives. The detail 
about model ([T]) can be found in Ref. [l5]. 



B. Spatial model 

Next we intend to add the spatial parts. Upto the first approximation, the dispersal of 
individuals can be taken random, so that Fick' law holds. This gives the flux terms as 



(3) 



where V 2 (V 2 = ^ + J^-) is the Laplacian operator in Cartesian coordinates. D s , D h 
and D r are the diffusion coefficients of the susceptible, infective, and recovered, respectively. 
Incorporating spatial terms into Eq. (CQ), it becomes 



— = A-dS - XSI + D S V 2 S, 
dt 



(4a) 

(4b) 
(4c) 



81 

— = \SI-{d + 7 )J - h(I) + AV 2 /, 
8R 

— = 7 J + h(I) -dR + D r V 2 R. 

Generally, we concern on the susceptible and infective individuals. Moreover the Eqs. (j4a| 
and (14b [1 are independent of the Eq. ( )4c|) whose dynamic behavior is trivial when /(to) — 
for some to > 0. So it suffices to consider the Eqs. ( |5al) and (|5bl) with I > 0. Thus, we 
restrict our attention to the following reduced spatial model 



<99 

— = A - dS - XSI + D S V 2 S, 
ot 

81 

— = \SI-(d + 1 )I-r + DiV 2 I. 



(5a) 



(5b) 



It is assumed that all the parameters are positive constants from the biological point of view. 



III. THEORETICAL ANALYSIS OF SPATIAL PATTERNS AND RESULTS 



To study the mechanism of the formation of spatial patterns, firstly, we analyze the sta- 
bility criterion of the local system. This can be obtained from the Ref. 15| . The system ([5]) 
has two positive equilibrium points if Rq > 1 and < H < (^/Rq — l) 2 , where Rq — XA 



l ° - d( rf+7 ) 
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and H = d ^ d ^ ■ The two positive equilibrium points are E\ = (S\,Ii) and E 2 = (5*2,12), 
where 

h = ^(Ro -l-H- ^{Ro-l-Hy-lH), 

51 = A/(d + \h), 

h = ^(Ro -1-H+ y/(Ro - 1 - Hf - AH), 

5 2 = A/(d + XI 2 ). 

Diffusion is often considered a stabilizing process, yet it is the diffusion-induced instability 
in a homogenous steady state that results in the formation of spatial patterns in a reaction- 
diffusion system [1]. The stability of any system is expressed by the eigenvalues of the 
system's Jacobian Matrix. The stability of the homogenous steady state requires that the 
eigenvalues have negative real parts. To ensure this negative sign, the trace of the Jacobian 
matrix must be less than zero at steady state if the determinant is greater than zero. 

The Jacobian matrix of system ([I]) at (S 2 ,I 2 ) is 

h=\~ d ~ Xh ~ XS2 |. (6) 
A/2 AS2 — d — 7 



From the Ref . 43j , we easily know that there are the Turing space in the system ([5]) at point 



E 2 , but at point E\ there is no Turing space. 



A. Stability of the positive equilibrium point in the spatial model 

In contrast to the local model, we employ the spatial model in a two-dimensional (2D) 
domain, so that the steady-state solutions are 2D functions. Let us now discuss the stability 
of the positive equilibrium point with respect to perturbations. Turing proves that it is 
possible for a homogeneous attracting equilibrium to lose stability due to the interaction 
of diffusion process. To check under what conditions these Turing instabilities occur in the 
model (jSJ), we test how perturbation of a homogeneous steady-state solution behaves in 
the long-term limit. Here we choose perturbation functions consisting of the following 2D 
Fourier modes 
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s = exp((k x x + k y y)i + 5 k t), (7a) 
i = exp((k x x + k y y)i + 6j,t). (7b) 

Since we will work with the linearized form of Eq. ((Sj) and the Fourier modes are orthog- 
onal, it is sufficient to analyze the long-term behavior of an arbitrary Fourier mode. 

After substituting S = S 2 + s and I = I 2 + i in Eq. (jSJ), we linearize the diffusion terms 
of the equations via a Taylor-expansion about the positive equilibrium point £^(5^, I 2 ) and 
obtain the characteristic equation 

= , (8) 





with 

J21 

here jn = -d - XI 2 , jn = —XS 2 , J21 = A/ 2 , and j'22 = AS*2 — d — 7. k 2 = k 2 x + k 2 y and k 
represents the wave numbers. 

To find Turing instabilities we must focus on the stability properties of the attracting 
positive equilibrium point E 2 (S2,l2)- The loss of stability occurs if at least one of the 
eigenvalues of the matrix J sp — 5^1 crosses the imaginary axis. From the Eqs. (jSJ) and (Q, 
we can obtain the characteristic equation 

det(J sp - 6 k I) =61- tr(J sp )4 + det(J sp ) = 0, (10) 

where tr(J sp ) = tr(J 2 ) - (D s + D t )k 2 and det(J sp ) = det(J 2 ) - fc 2 (j u A + j 22 D s ) + k A D s D % . 
Taking tr(J 2 ) > tr( J sp ) into account, we can obtain that for saddles and attractors (both 
with respect to the non-spatial model) a change of stability coincides with a change of the 
sign of det(J sp ). 

Doing some calculations we find that a change of the sign of det(J sp ) occurs when k 2 
takes the critical values 



, 2 = ./i.M • nil I- - y/(jnA + J22^s) 2 - 4P a Adet(J 2 ) . , 

2D s D l ' 1 aj 

k 2 = JuDj+jM + VUuDi + J22A) 2 - 4AAdet(J 2 ) ( 
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FIG. 1: (Color online) This graphes illustrate the eigenvalues of the spatial model ([5]) at positive 
equilibrium point ~E^i{S% % I2), and the loss of stability occurs relation to the limit range wave 
numbers, (a) For diffusion-driven instability arise both tr(J sp ) and det(J ap ) must be negative for 
some range of k 2 ; (b) The real parts of the eigenvalues of matrix (JH) at positive equilibrium (S2, ^2)- 
Model parameters used here are: A = 3, d = 0.3, A = 0.35, r = 0.5, 7 = 0.8, D s = 0.02, and 
Di = 0.0005. 



In particular, we have 

det(J sp ) < kl < k 2 < k\. (12) 

If both k 2 _ and k 2 + exist and have positive values, they limit the range of instability for a 
local stable equilibrium. We refer to this range as the Turing Space (or Turing Region, see 

Fig. CD . 

In Fig. CD, the real parts of the eigenvalues of the spatial model (JSJ) at positive equilibrium 
point E 2 (S2, h) are plotted. From the Eqs. ( J7a|) and ( l7bl) . we know that the parameter 
5k can either be a real number or a complex number. If it is a real number, the spatial 
patterns will emerge and be stable over time and otherwise the spatial patterns will die 
away quickly. In both case, the sign of the real parts of 5k (written as Re(5k)) is crucially 
important to determine whether the patterns will grow or not. In particular if Re(5k) > 0, 
spatial patterns will grow in the linearized system because le^l > 1, but if Re(5k) < the 
perturbation decays because \e Sk \ < 1 and the system returns to the homogeneous steady 
state. Further details concerning linear stability analysis can be found in Ref. 44j. The 
Fig. [1] presents the typical situation of Turing instability With respect to homogenous 
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FIG. 2: (Color online) The graphes illustrate Turing space versus the parameter A and r, re- 
spectively. In the Turing space fc_ and k + curves are shown respectively, (left) The Bifurcation 
diagrams show the formation of the two stationary solutions of Eq. © with fixed r and varying A, 
and values of the parameters are A = 3, d = 0.3, r = 0.5, 7 = 0.8, D s = 0.02, and Di = 0.0005; In 
the lower diagram, the green dashed line represents the loss of stability for the positive equilibrium. 
The red line represents the stable equilibrium, (right) The values of the parameters are A = 3, 
d = 0.3, A = 0.35, 7 = 0.8, D s = 0.02, and A = 0.0005. 

perturbations, E 2 (S 2 ,I 2 ) is stable at first, but when k 2 increases, one eigenvalue changes its 
sign (when k 2 arrives at k 2 _.), the instability occurs. The instability exist until k 2 reaches 
k 2 ^. When k 2 is over k 2 +) (S 2 ,I 2 ) returns stability again. Thus the Turing space is bounded 
between k 2 _ and k\. 

The change of the bounds A;_ and k + with respect to the variation of the A and r are 
illustrated in the Fig. El respectively. The typical feature of Turing space in the model (jSJ) can 
be observed in Fig. El The Turing space is limited by two different bounds. From Fig. El^ left- 
top) we can see that and k + converge in one point A which corresponds to the critical 
value, A c . Beyond right bound (line a), the E 2 (S 2 ,l2) exists and is stable. The left bound 
(line b) of the Turing space shows an "open end", which corresponds to the saddle-node in 
the bifurcation plot (see Fig. El(left-bottom)) for the model (CE]) and the equilibrium point 
E 2 {S 2 il2) does not exist under this bound. This figure shows the solutions of equilibrium 
points 1 1 and I 2 , where the solid curves represent attractors, dashed curve represents the 
repellers and saddles, the dotted line a represents the periodic points. This diagram explains 
the E 2 (S 2 ,I 2 ) changes from repeller to attractor, and an unstable orbit of periodical points 
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FIG. 3: (Color online) Basic dispersion relation giving the growth rate Re(5k) as a function of the 
wavenumber k 2 . The mode become marginal at the point (H, A = A c ) a finite-/c 2 (Turing) mode. 
The parameter values are A = 3, d = 0.3, r = 0.5, 7 = 0.8, D s = 0.02, and D, L = 0.0005. 

emerges. From the Fig. fright), the curves indicates that and k + converge in one point 
(A). Below that bound, the E 2 (S 2 ,l2) exists and is stable. The right bound of the Turing 
space also shows an "open end" . 

Comparing to the two graphes in Fig. [2], one can obtain that the parameters A and r 
have a similar dynamical behavior in the system for the Turing-bifurcation, but their effects 
are opposite. We use the parameter A as the Turing-bifurcation parameter in present paper. 
Fig. [3] shows growth rate curves of the spatial patterns, where at bifurcation (curve b), 
A = A c ~ 0.547 (The threshold A c can be derived analytically, see the Appendix. The 
analytical and numerical values of A c are approximately equal.), from spatially uniform to 
spatially heterogenous the critical wave number (point H) is k c = ^Jk 2 c + k 2 c . The curves 
a and c correspond to the parameter A = 0.35 below the A c , and A = 0.8 above the A c , 
respectively. The spatial patterns are generated when A passes through the critical Turing- 
bifurcation point A c . And for A < A c there is a finite range of unstable wave numbers which 
grow exponentially with time, 0(exp(Skt)), where 5^ > for a finite range of k. 

The stable characteristics of E^OS^,-^) can be changed by the variation of parameter A: 
A sufficiently high increase of A will turns E 2 (S 2 ,l2) into an attractor. When changing its 
characteristics, E 2 (S 2 ,l2) traverses a subcritical Hopf bifurcation and an unstable periodic 
orbit emerges (the dotted line a in Fig. EJ^left)). Surprisingly, the latter is not necessarily 
true, if effect of diffusion comes to play. 
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B. Spatial patterns of the spatial model 



The numerical simulations are performed in this section for the spatial model (jSJ) in two 
dimensions. During the simulation, the periodic boundary conditions are used and part 



of the parameter values can be determined following Ref. [15( (see the Fig. [I] and [2]). We 
assume that the homogeneous ^(S^, I2) distributions are in uniform states for each start of 
the simulation. To induce the dynamics that may lead to the formation of spatial patterns, 
we perturb the /-distribution by small random values. 




FIG. 4: (Color online) Snapshots of contour pictures of the time evolution of I(x,y,t) at different 
instants, (a)-(c) Numerical results in 100 x 100 sites. The parameter values are A = 3, d = 0.3, 
r = 0.5, 7 = 0.8, A = 0.65, D s = 0.02, D { = 0.0005, and Ax = Ay = 0.05. (A) iteration; (B) 
5000 iterations; (C) 30000 iterations; (D) 40000 iterations. [Additional movie formats available 
from the author] 

We study the spatial model (jSJ) by performing stable analysis of the uniform solutions and 
by integrating Eqs. ( |5al) and ( I5bl) numerically at different values of A on a grid of 100 x 100 
sites by a simple Euler method with a time step of At = 0.01. The results for the infected 
are summarized below in two dimensions. The model has a uniform free-disease state (no 
infected) for all constant values of A, represented by the solution S = A/d, 1 = 0. The free- 
disease state is stable when A < A c /. Here X c > is a critical value (or threshold ) corresponded 
by the dotted line b in Fig. [2] (left). Above A c / two new states appear, shown in Fig. [2] 
(left) as line I\ and I2. The state I2 represents a uniformly distributed population with 
infected density monotonically increasing with A. It is instable only for relative values of A, 



10 




FIG. 5: (Color online) Snapshots of contour pictures of the time evolution of I(x,y,t) at different 
instants. The parameters values are the same as Fig. HI (A) iteration; (B) 5000 iterations; (C) 
40000 iterations; (D) 42000 iterations. [Additional movie formats available from the author] 




FIG. 6: (Color online) Snapshots of contour pictures of the time evolution of I(x,y,t) at different 
instants. The parameters values are the same as Fig. U] but A = 0.5. (A) iteration; (B) 1500 
iterations; (C) 8500 iterations; (D) 9400 iterations. [Additional movie formats available from the 
author] 

A c ' < A < A c and regains stability when A > A c , where the infected density is high. The 
types of spatial patterns are depended on the range of parameter A as in Figs. HJ El and 
We have made movies of Figs. HJ [5] and El respectively, as supplementary materials. 

We test several different initial states within the linear regime and the nonlinear regime 
respectively. Figs. H] and [5] show that stationary stripe and spot patterns emerge mixtedly 
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in the distribution of the infected population density, where the A is above A c in the linear 
regime. The initial state of Fig. H] is the random perturbation of the stable uniform infected 
state. The initial state of Fig. consists of a few spots (100 scattered spots). Values of the 
parameters are the same in both two figures. In the linear regime, the result shows that 
the stripes and spots which describe asymptotic patterns for the spatial model §5§ converge 
at long times. Different initial states may lead to the same type of asymptotic patterns, 
but the transient behaviors are obviously different (compare the Fig. IH(b) with Fig. [5](b)). 
Unfortunately, the linear predictions are not accurate in the nonlinear regime. 

Fig. [6] shows the snapshots of spatial patterns when the A is between X c i and A c . One can 
see that the spatial patterns resemble in Fig. @] and Fig. E] at the beginning phase, but differ 
in the middle and last phases. In Fig. [HI the spatial patterns appear spotted, holed and 
labyrinthine states in the middle phases, and the spatial patterns appear uniform spatial 
states in the last phase. 

To explain spatial patterns arising from the spatial model, here we present some observa- 
tions of the spatial and temporal dynamics of dengue hemorrhagic fever epidemics. Dengue 
fever (DF) is an old disease that became distributed worldwide in the tropics during the 
18th and 19th centuries. Fig{7| shows spotted and labyrinth- like spatial patterns of DHF 



from the field observations 



451 ] . By comparing the Figs. H] and [5] with Fig. [7J our results 



simply capture some key features of the complex variation and explain the observation in 
spatial structure to most vertebrate species, including humans. In Figs. 0] and O the steady 
spatial patterns indicate the persistence of the epidemic in the space. This result well agrees 
with the field observation. More examples of the spatial patterns of epidemic, such as HIV, 



poliovirus, one can find in Refs. 46|, |47|. In the light of recent work of emphasizing the exis- 



tence of 'small world' networks in human population, our results are also consistent with M. 
Boots and A. Sasaki's conclusion that if the world is getting 'small'-as populations become 
more connected-disease may evolve higher virulence 



IV. DISCUSSION 



In our paper, the labyrinthine patterns are found in the spatial epidemic model driven by 
the diffusion. The spatial epidemic model comes from the classical non-spatial SIR model 
which assumes that the epidemic time scale is so tiny related to the demographic time scale 
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FIG. 7: (Color online) Field observations of DHF spatial patters. Reprinted form Trends In 



Microbiology, Ref. 



451 ] . Copyright (2002), with permission from Elsevier. Aedes aegypti distribution 



in the Americas: 1930s, 1970 and 2001. 



that demographic effects may be ignored. But here in the spatial model we take the births 
and deaths into account. The spatial diffusive epidemic model is more realistic than the 
classical model. For instance, the history of bubonic plague describes the movement of 
the disease from place to place carried by rats. The course of an infection usually cannot 
be modeled accurately without some attention to its spatial spread. To model this would 
require partial differential equations (PDE), possibly leading to descriptions of population 
waves analogous to disease waves which have often been observed. Here our spatial model 
is established from a basic dynamical 'landscape' rather than other perturbations, including 
environmental stochastic variations. From the analysis of the Turing space and numerical 
simulations one can see that the attracting positive equilibrium will occur instability driven 
by the diffusion and the instability leads to the labyrinthine patterns within the Turing 
space. This may explain the prevalence of disease in large-scale geophysics. The positive 
equilibrium is stable in the non-spatial models, but it may lose its stability with respect 
to perturbations of certain wave numbers and converge to heterogeneous distributions of 
populations. It is interesting that we have not observed the isolated spots patterns in the 
spatial epidemic model (jSJ). 

The model (jHJ) is introduced in a general form so that it has broad applications to a range 
of interacting populations. For example, it can be applied to diseases such as measles, AIDS, 
flu, etc. Our paper focuses on the deterministic reaction-diffusion equations. However, recent 



13 



study shows noise plays an important role on the epidemic model [48l . |49j, which indicate 
that the noise induces sustained oscillations and coherence resonance in the SIR model. 
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APPENDIX A: APPENDIXES 

Considering the dispersive relation associated with Eqs. ({Taj) and fl7b|) . the functions of 
8k for the spatial model are defined by the characteristic Eq. (flUj) . Now Re(5k) predict the 
unstable wave modes from Eq. fflOl) . One can estimate the most unstable wave number and 
the critical value of the bifurcation parameter by noticing that at the onset of the instability 
5k{k c ) = 0. Thus the constant term in Eq. ( TlDl) must be zero at k c . In the case of the spatial 
model this condition is a second order equation on k 2 , i.e., 

D s D % k\ - k 2 c ((-d - XI 2 )D t + (XS 2 -d- j)D s ) 

+(-d - XI 2 )(XS 2 - d - 7 ) + X 2 S 2 I 2 = 0. (Al) 

And as a result the most unstable wave number is given by (~ d ~ XI ^) D ^+( x ^~ d ~~i) Ds _ ^he 
critical Turing-bifurcation parameter value, which corresponds to the onset of the instability 
is defined by Eq. ( 1A1I) . In the spatial model A is the bifurcation parameter adjusting the 
distance to the onset of the instability. The discriminant of Eq. (IA1I) equals zero for A c and 
an instability takes place for A < A c . Then, we have 

A (D S S 2 + Dih + 2^D s D i S 2 I 2 ){D sl + D s d - Djd) 

(Aia " D s S 2 f 

(A2) 

where S2 = S 2 (X C ) and I 2 = I 2 (X C ). We can calculate A c from the Eq. (1A2I) by the computer. 
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